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Abstract 

We present the implementation of a solution scheme for fluid-structure interaction prob- 
lems via the finite element software library deal. II. The solution scheme is an immersed 
finite element method in which two independent discretizations are used for the fluid and 
immersed deformable body. In this type of formulation the support of the equations of mo- 
tion of the fluid is extended to cover the union of the solid and fluid domains. The equations 
of motion over the extended solution domain govern the flow of a fluid under the action of 
a body force field. This body force field informs the fluid of the presence of the immersed 
solid. The velocity field of the immersed solid is the restriction over the immersed domain 
of the velocity field in the extended equations of motion. The focus of this paper is to show 
how the determination of the motion of the immersed domain is carried out in practice. We 
show that our implementation is general, that is, it is not dependent on a specific choice 
of the finite element spaces over the immersed solid and the extended fluid domains. We 
present some preliminary results concerning the accuracy of the proposed method. 

Keywords: Fluid Structure Interaction; Immersed Boundary Methods; Immersed Finite 
Element Method; Finite Element Immersed Boundary Method 



1 Introduction 

Heltai and Costanzo (2012) have recently discussed a fully variational formulation for an im- 
mersed method to the solution of fluid-structure interaction (FSI) problems. Immersed methods, 
which deal with the motion of bodies immersed in fluids, allow one to choose the discretization 
for the fluid and solid domains independently from each other. As such, they stand in contrast to 
established methods like the arbitrary Lagrangian-Eulerian (ALE) ones (see, e.g., Hughes et al., 
1981), where the topologies of the solution grids for the fluid and the solid are constrained. 
Immersed methods have three main features: 
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1. The support of the equations of motion of the fluid is extended to the union of the physical 
fluid and solid domains. 

2. The equations of motion of the fluid have terms that, from a continuum mechanics view- 
point, are body forces "informing" the fluid of its interaction with the solid. 

3. The velocity field of the immersed solid is identified with the restriction to the solid domain 
of the velocity field in the equations of motion of the fluid. 

A taxonomy of immersed methods can be based on how these three elements are treated theo- 
retically and/or are implemented practically (see the discussion in Heltai and Costanzo, 2012). 
Here we employ the approach proposed in Heltai and Costanzo (2012) in which the entire so- 
lution scheme is developed within the general framework of the finite element method. Most 
importantly, the restriction mentioned at point 3 above is done via a fully variational approach. 
As such, the approach demonstrated herein stands in contrast to what is used in the immersed 
boundary methods stemming from the approach of Peskin and his co-workers (see, e.g., Griffith, 
2012; Griffith and Luo, 2012; Peskin, 1977, 2002) or the finite element extension of Peskin's 
approach due to Liu and co-workers (see, e.g., Liu et al., 2007; Wang and Liu, 2004; Zhang 
et al., 2004), which is based on the implementation of the reproducing kernel particle method. 
As explained in detail in Heltai and Costanzo (2012), the method demonstrated herein stems 
from the approach by Boffi and Gastaldi (2003), Heltai (2006), and Boffi et al. (2008), and Heltai 
(2008). 

In Section 2 we review the problem's governing equations. In Section 2.3 we present the 
variational reformulate of the governing equations and we will present their discrete counterparts 
in Section 3. The content of Sections 2 and 3 follows closely the exposition in Heltai and Costanzo 
(2012) and is reported here for completeness. In Section 4 we provide details about the code we 
have developed and instructions for compilation, execution, and generation of documentation. 
The entire code is based on the open source deal. II library (see Bangerth et al., 2006, 2007). 
We conclude the article with Section 5, where we present some numerical results. 

2 Problem Formulation 

2.1 Basic notation and governing equations 

B t in Fig. 1 represents the configuration of a regular body at time t. B t is a (possibly multiply 
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Figure 1: Current configuration Bt of a body SS immersed in a fluid occupying the domain f2. 

connected) proper subset of a fixed control volume $7. The domain O \ Bt is filled by a fluid 
and we refer to Bt as the immersed body. dO, and dBt, with outer unit normals m and n, 
respectively, are the boundaries of Q and Bt. We denote by B the reference configuration of the 
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immersed body. We denote the position of points of SB in B by s, whereas we denote the position 
at time t of a generic point P G f2 by xp(t). A motion of S3 is a diffeomorphism £ : B — >■ 
a; = £(s,t), with s g 5, i £ fl, and t € [0,T), with T a positive real number. 

The function p(x, t) describes the mass density in the entire domain fl. The function p can 
be discontinuous across dB t . The local form of the balance of mass requires that, Vi G (0,T), 

p + pdivw = 0, aj€n\(8ftU0Bt), (1) 

where u(x,t) = d£(s,t)/dt\ g=( >- 1 , t -. is the velocity field, a dot over a quantity denotes the 
material time derivative of that quantity, and where 'div' represents the divergence operator 
with respect to x. 

The local form of the momentum balance laws require that, Vt G (0,T), T = T T (the 
superscript T denotes the transpose) and 

divT + pb = pii, xen\{dQUdB t ), (2) 

where J(x,t) is the Cauchy stress and b(x,t) is the external force density per unit mass acting 
on the system. 

In addition to Eqs. (1) and (2), we demand that the velocity field be continuous (correspond- 
ing to a no slip condition between solid and fluid) and that the jump condition of the balance 
of linear momentum be satisfied across dB t : 

u(x + , t) = u(x~ , t) and T(x + , t)n = T(x~, t)n, x G dBt, (3) 

where the superscripts — and + denote limits from within and without Bt , respectively. 

We denote by dfl, d and 8Qn the subsets of dil where Dirichlet and Neumann boundary data 
are prescribed, respectively. The domains dVL^ and dVLjv are such that 

dQ = dQ D UdQ N and dSl D ndn N = (4) 

We denote by u g (x,t), with x G 8Qd, and by T g (x,t), with x G dtljv, the prescribed values of 
velocity (Dirichlet data) and traction (Neumann data), respectively, i.e., 

u(x,t) = u g (x,t), for x G 8£Id, and T(x, t)m(x, t) = T g (x, t), for x G <9f2jv> (5) 

where the subscript g stands for 'given'. 

2.2 Constitutive behavior 

Constitutive response of the fluid. We assume that the fluid is linear viscous and incom- 
pressible with uniform mass density p. Denoting by Tf the constitutive response function of the 
Cauchy stress of the fluid, we have (see, e.g., Gurtin et al., 2010) 

t f = -pl + 2 M D, D = i(L + L T ), (6) 

where p is the pressure of the fluid, I is the identity tensor, p > is a given viscosity coefficient, 
and L = gradu, and where a "hat" (T) is used to distinguish the constitutive response function 
for T from T itself. For convenience, we denote by T^ the viscous component of Tf, i.e., 

T f / = 2 / uD = ^(L + L T ). (7) 

Due to incompressibility, the balance of mass reduces to 

divw = fovxen\B t . (8) 

Under these conditions, p is a Lagrange multiplier that allows to enforce Eq. (8). 
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Constitutive response of the solid. The immersed body is taken to be incompressible and 
viscoelastic of differential type: 

ts = -p\ + f e B + f v B , (9) 

where T| and Tg denote the elastic and viscous parts of T s , respectively, and p is the Lagrange 
multiplier needed to enforce incompressibility. The viscous part of the behavior is assumed to 
be of the same type as that of the fluid, that is, 

t^ = 2/xD = /x(L+L T ), (10) 

where /x is the same constant viscosity coefficient of the fluid. We assumed that Tg derived from 
a strain energy potential. To be precise, let the first Piola-Kirchhoff stress tensor be P. This 
tensor is related to T as follows (see, e.g., Gurtin et al., 2010): 

P = JTF- T , (11) 

where J = det F, and the tensor F, called the deformation gradient, is defined as 

F-« (12) 

OS 

Letting P| = JT^F denote the constitutive response function for the elastic part of the first 
Piola-Kirchhoff stress tensor, as is typical in elasticity, we assume the existence of a function 
W" s e (F) such that 

dW e (F) 

Ps = (13) 

where Wg is the density of the elastic strain energy of the solid per unit volume. Invariance 
under changes of observer demands that W B be a function of an objective strain measure such 
as C = F T F. If the solid is isotropic, W B must be a function of the principal invariants of C. 

2.3 Reformulation of the governing equations 

We now reformulate the governing equations in variational form. The motion of the solid will 
be described via the displacement field, denoted by w and defined as 

w(s,t) := C(s,t) - s, s£B. (14) 

The displacement gradient relative to the position in B is denoted by H: 

H:=£ * H-F-.. (!5) 

Equation (14) implies 

w(s,t) = u(x,t)\ x=c{sty (16) 

The principal unknowns of our fluid-structure interaction problem are then the fields 

u(x,t), p(x,i), and w(s,t), with x G 0, s G B, and t G [0, T). (17) 

The functional spaces for the problem are 

uer = Hjyinf := {u G L\Q) d | V x u G L 2 (Q) dxd , u\ d n D = (18) 

pe £ := L 2 (fl), (19) 

w G & = H 1 (B) d := <j™ G L 2 (B) d \ V s w G L 2 (B) dxd ^, (20) 

where W x and V s denote the gradient operators relative to x and s, respectively. Also, referring 
to Eq. (18), the function space for the test functions for the velocity field is taken to be as 
follows: 

% = Hl(Sl) d := [v e L 2 (n) d \\/ X v e L 2 {n) dxd ,v\ 9QD =o}. (21) 
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2.4 Variational restatement of the governing equations 

When the solid is incompressible, the mass density of both the fluid and the solid are constant 
so that p = (almost) everywhere in 0,. Then, referring to Eqs. (5), Eqs. (18)-(20), and 
the constitutive response functions of both the fluid and the solid, the governing equations 
introduced so far can be expressed in weak form as follows: 

/ p(u — b) ■ v dv + / Tf • V x v dv 
Jn Jn 

+ / (t s -T{)-V x vdv- r g -vda = Mve% (22) 
JB t Jan N 

and 

qdivudv = Vg G B. (23) 

Jn 

A crucial aspect of our approach is the enforcement of Eq. (16). We enforce this relation weakly 
as follows: 

<Z> B j \w(s, t) - u(x, t) \ x=c{S:t) ] • y(s) dV = Vy G &, (24) 

J B 

where dV is an infinitesimal volume element of B, and where $b is a constant with dimensions 
of mass over time divided by length cubed, i.e., dimensions such that, in 3D, the volume integral 
of the quantity has the same dimensions as a force. We observe that, since we have 

assumed that the viscous part of the stress response of the solid is the same as that of the fluid 
(Heltai and Costanzo, 2012 discuss the most general of cases in which the immersed body and 
the surrounding fluid can have different constitutive response functions), the term (T s — Tf) in 
Eq. (22) is equal to the elastic response of the solid T|. 

Our numerical approximation scheme for Eqs. (22)-(24) is based on the use of two inde- 
pendent triangulations, namely, one of O and one of B. The fields u and p, as well as their 
corresponding test functions, will be expressed via finite element spaces supported by the trian- 
gulation of Q. By contrast, the field w will be expressed via a finite element space supported by 
the triangulation of B. Because of this, any term in Eq. (22) defined over B t is now rewritten 
as an integral over B: 

/ p{ii — b) ■ v dv — / pdivvdv+ / Tf-V x vdv— / r g -vda 
Jn Jn Jn Jan N 

+ ! P e s F T (s,t)-V x v( X )\ x=c{st) dV = Vve%. (25) 

We now define the operators we will use in our finite element formulation. In these definitions, 
we will use the following notation: 

v^^v, (26) 

in which, given a vector space V and its dual V* , ip and 4> are elements of the vector spaces V* 
and V, respectively, and where v ,(»,»^ v identifies the duality product between V* and V. For 
convenience, we also introduce the following shorthand notation 

t>] = p[V x u(x, t) + (V x u(x, t)) T ] , (27) 
F[w] = I + V s w(s, t), (28) 
dWHF 



dF 



(29) 

F=F[w] 



Finally, to help identify the domain and range of these operators, we establish the following 
convention. We will use the numbers 1, 2, and 3 to identify the spaces Y, J3, and respectively. 
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We will use the Greek letter a, (3, and 7 to identify the spaces Y*, J3* , and <3f* , respectively. 
Then, a Greek letter followed by a number will identify an operator whose domain is the space 
corresponding to the number, and whose co-domain is in the space corresponding to the Greek 
letter. For example, the notations 

£ a2 and £ a2 p (30) 

will identify a map (£ a 2) from £! into Y* and the action of this map (£ a 2P £ Y*) on the field 
p £ J, respectively. If an operators has only one subscript, that subscript identifies the space 
containing the range of the operator. With this in mind, let 

r 



Mai 


■Y- 


» Y* 


Nal{u) 


: Y - 


■»• Y* 


Vol 


: Y - 


■»• Y* 


% 


: Y - 





(A1 a itt,v) r := pu-vdv Vu£Y,\/v £ Y , (31) 

(A/"ai(w)w,v) r := / ^(V^Ju • -u dv \/u,w £Y,\/v£ Yo, (32) 

r ,(Paiw,v> r := / f^M-V^du V M ef,VDer , (33) 
Jn 

^{B pi u,q)^:=- qdivudv \/q£^,\/u£Y. (34) 

The operators defined in Eqs. (31)-(33) are found in traditional variational formulations of the 
Navier-Stokes equations and will be referred to as the Navier-Stokes component of the problem. 
As typical of other immersed methods, these operators have their support in ft as a whole. 

We now define the operator in our formulation that has its support over B but does not 
contain prescribed body forces or boundary terms. 

A a (w,h) G Y*, Vto, h G ^,Vw Ef,V« Y 

r.(Aa(w,h),v) y := / [P S >]F>] • V x v(x)] x=a+h(s , t) dV. 

J B 

We now define operators with support in B that express the coupling of the velocity fields 
defined over Vl and over B. Specifically, we have 

M y3 ■ & -> &*, Vw,y£ 

. . f (36) 

&*{M 7 3W,y)g, := $ B / w-y(s)dV, 

Mjiiw) : Y -> V«£ Y,Vw,y£ & , 

y . <M 7 i («;)«, y> y :=$s / u(a^)| s=s+it)(s t) -y(s)dy, 

Finally, we define the operators that express the action of prescribed body and surface forces. 

F a g y\ vb g H-\n),vr g g H~^(dn N ), g r 

. . /■ f (38) 

y*(T a ,v) y := pbvdv + r g -vda 
Jn Jan N 

G a {w) G r*, Vu> G ^,V6 G i7 _1 (n),Vi; G Y 

r .(g a (w),«> r := / (p S0 (s)- P J[w])b-v{x)\ x=g+w{st) dv. 

J B 

In the definition of the operator A a in Eq. (35), the motion of the immersed solid plays a 
double role in that it affects the elastic response of the solid (through w) as well as the map 
(through h) functioning as a change of variables of integration. As discussed in Heltai and 
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Costanzo (2012), it is important to separate these two roles and view A a as the composition of 
a change of variable operator and a Lagrangian elastic operator. To do so, we write 

S ai (h) : ->• r*, Vy* G J^,Wi G ^,V^ G r 

r(Sc«(h)y*,v)y := .^*{y* Mx)\ x=s+h{s) ) ^ 1 U " 

Ay{w) G JT y *, Vto G r,Vy G 

(41) 

( /.) [ill I 'J J \ / H' lip . \ n i i \ 

■y££> * 



..<AM,y>^ := / P|H-V s ydV. 



Once the operators S ai (K) and Ay(w?) are defined, one can prove the following theorem (see 
Heltai and Costanzo, 2012): 

Theorem 1 (Eulerian and Lagrangian elastic stiffness operators of the immersed domain). With 
reference to the definitions in Eqs. (35), (40), and (41), we have 

A a (w,h) =S on (h)A y {w) and S aj (h) = M^i(h)M~^, (42) 

where S ai (h)Aj(w) and A4^i(h)A4~^ indicate the composition of the operators S ai (h) and 
Ay (to) and of the operators Ai^j^h) and M.~l , respectively. 

The operators defined above allow us to formally restate the overall problem described by 
Eqs. (25), (23), and (24) as follows: 

Problem 1 (Dual formulation). Given initial conditions Uoff and Wq G & , for all t G (0, T) 
find u(x,t) G y, p(x,t) G £?, and w(s,t) G such that 

M al u' + M al (u)u + V al u + (B jS i) T p + S ai {w)A 1 {w) =F a + Qa{w), (43) 

Bpiu = 0, (44) 
Mysw' - M 7 i(w)u = 0, (45) 

where u'(x,t) = du(x,t)/dt and w'(s,t) = dw(s,t)/dt. 

Problem 1 can be formally presented in terms of the Hilbert space 3? := V x =S x and 
3?q := Yo x =S x with inner product given by the sum of the inner products of the generating 
spaces. Defining 3? 3 £ := [u,p, w] T and 3fo 3 i[) := [v, q, y] T , then Problem 1 can be compactly 
stated as 

Problem 2 (Grouped dual formulation). Given an initial condition £o £ for all t G (0, T) 
find f (t) G 3?, such that 

{T(t,£,?),i/>) = 0, V^GiTo, (46) 
where the full expression of T : 2f >->■ Jq* is defined as in Problem 1 . 

The energy estimates concerning the above abstract formulation has been discussed in Heltai 
and Costanzo (2012) where it is shown that stability is obtained under the same assumptions 
that yield stability for the Navier-Stokes problems. 



3 Discretization 

3.1 Spatial discretization by finite elements 

The fluid domain is discretized into the triangulation VL^ and the immersed body into the tran- 
gulation Bh ■ each of these triangulations consists of (closed) cells K (triangles or quadrilaterals 
in 2D, and tetrahedra or hexahedra in 3D) such that: 
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1. n = U{K G Q h }, and B = U{K G B h }; 

2. Any two cells K, K' only intersect in common faces, edges, or vertices; 

3. The decomposition matches the decomposition dft = dQp U 8Qn ■ 

On Qh and Bh, we define the finite dimensional subspaces % C Y , £2h C =S, and ^ C ^ as 
follows: 

*h := {u h G V | G * € = spaa{«* ft }S ( 47 ) 

<2 ft := | G P Q (if), Ken h } = span{^}£ (48) 

% := {w h e&\ w h \ K G Py(if), i^GB^ spaa{yU£i> (49) 

where Vv(K), Vq(K) and Vy{K) are polynomial spaces of degree ry, tq and ry respectively on 
the cells if, and iVy, ^Vq and Ny are the dimensions of each finite dimensional space. The pair 

and J2h are chosen so that the inf-sup condition for the well-posedness of the Navier-Stokes 
problem (see, e.g., Brezzi and Fortin, 1991) is satisfied. 

The discrete version of Problem 1 is now presented using a matrix notation. An element of 
a discrete space, say Uh G is represented by a column vector of time dependent coefficients 
u h(t)i 3 = 1) • • • j such that Uh(x, t) = ^ ul(t)vi(x) , where v J h is the j th base element of 
Yh- We use the notation M a ±Uh to represent the multiplication of the column vector Uh by the 
matrix whose elements are 

MSi-=AM a iv{,v\) y , (50) 

where the operator in angle brackets is the one defined earlier. A similar notation is adopted for 
all other previously defined operators. With this notation, the duality products in the discrete 
spaces are indicated by simple scalar products in ~K N (N depending on the dimension of the 
system at hand). Hence, using the matrix M a i, we can write 

y,(M a iu h ,v h ) y = v h -M al u h , (51) 

where the dot-product on the right hand side is the scalar product in R Ny . 

Chosen fi/, and Bh along with £?h, and &h, Problem 1 is reformulated as follows: 

Problem 3. Given u G %, w G for all t G (0,T), find u h (t) G %, Ph(t) G ^h, and 



Wh(t) G &h such that 

M al u' h + N a i(u h )u h + D al u h + {Bpxfph + S ai {w h )Ay(w h ) = F a + G a (w h ), (52) 

B pi u h = 0, (53) 

M j3 w' h -M 7l (w h )u h = 0, (54) 



where u' h (x,t) = J2[ u h(t)Y v h( x ) an d w 'h( 8 ^) = YH w h(t)]' Vhi 8 ) i and where the prime denotes 
ordinary differentiation with respect to time. 

In compact notation, Problem 3 can be casted as semi-discrete problems in the space 3? D 

Problem 4. Given an initial condition £o £ ^h, for all t G (0, T) find £h{t) G 3?h, such that 

F(t,£ h ,? h ) = 0, (55) 

where 

:=(T(t^ h ,e h ),^i), i = 0,...,N v + N Q + N Y , (56) 

and T has the same meaning as in Eq. (46), with ip h being the basis function for the spaces 
or 1¥h corresponding to the given value of i. 
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3.2 Coupling of the fluid and immersed domains 

The operators M a i, N a i(uh), D a i, Bpi, and F a in Problem 3 are common in variational formu- 
lations of the Navier-Stokes problem and were implemented in a standard fashion. The operator 
M 7 3 was also implemented in a standard fashion since it is the mass matrix for Less common 
are the operators that depend nonlinear ly on the motion of the immersed domain w. Thus, we 
now discuss the practical implementation of such operators. 

Let's consider, for example, the matrix M~x(w) contributing to the velocity coupling between 
the fluid and immersed domain: 

Mii 1 {w h )= {M yl {w h )vlyi) Wh =$ B f v j h (x)\ x=s+Wh{a>t) ■ y\(s) dV. (57) 

The above integral is computed by summing contributions from each cell K of Bh- Each of 
these contributions is a sum over the Nq quadrature points. We observe that the integrand 
y\X s ) is supported over the triangulation of Bh but the functions v ] h {x) (with x = s + Wh(s,t)) 
are supported over the triangulation Jl^. Therefore, the construction of operators like M l A[wh) 
draws information from two independent triangulations. In our code, we start by determining 



• : quadrature point of the solid cell 



fluid cell C 
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Figure 2: Cells denote as A-D represent a four-cell patch of the triangulation of the fluid domain. 
The cell denoted as "solid cell" represents a cell of the triangulation of the immersed solid domain 
that is contained in the union of cells A-D of the fluid domain. The filled dots represent the 
quadrature points of the quadrature rule adopted to carry out integration over the cells of the 
immersed domain. 

the position of the quadrature points of the immersed element, both relative to the reference 
unit element and relative to the global coordinate system adopted for the calculation, through 
the mappings: 

s K : K := [0, l] d ^ K G B h , (58) 
I + w h : K ^ solid cell. (59) 

These maps allow us to determine the global coordinates of the quadrature points. These 
coordinates are then passed to a search algorithm that identifies the cells in Qh that contain 
the points in question. In turn, this identification allows is to evaluate the functions vi . The 
overall operation is illustrated in Fig. 2 where we show a cell of Bh straddling four cells of f2/j 
denoted fluid cells A-D. The quadrature points over the solid cell are denoted by filled circles. 
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The contribution to the integral in Eq. (57) due to the solid cell is then computed by summing 
the partial contributions corresponding to each of the fluid cells intersecting the solid cell in 
question: 

K£B h JK 

~ E Y.<^ x )\ x =s K , q+Wh {s K , q ,t)-Vh^>K, q , (60) 

where sx,q is the image of q-ih. quadrature point under the mapping sk, and 0JK,q is the corre- 
sponding quadrature weight. The implementation of an efficient search algorithm responsible for 
identifying the fluid cells intersecting an individual solid cell is the only technically challenging 
part of the procedure. We use the built-in facilities of the deal . II library to perform this task. 
Once the fluid cells containing the quadrature points of a given solid cell are found, we determine 
the value of v J h at the quadrature points using the interpolation infrastructure inherent in the 
finite element representation of fields defined over Jl^. The deal. II C++ class we use for this 
implementation is the FEFieldFunction. 



3.3 Time discretization 

Equation (55) represents a system of nonlinear differential algebraic equations (DAE), which we 
solve using a Newton iteration. In the code accompanying this paper, the time derivative £' is 
approximated very simply via an implicit-Euler scheme: 

i'n = h-\tin-in-l), (61) 

where and £' n are the computed approximations to £(t n ) and £'(t n ), respectively, and the step 
size h = t n — t n -\ is kept constant throughout the computation. Although not second order 
accurate, this time stepping scheme is asymptotically stable. 

The application of the implicit-Euler scheme in Eq. (61) to the DAE system in Eq. (55) 
results in a nonlinear algebraic system to be solved at each step: 



G(&) :=F{t n ,^h~ l (i n -i n ^ =0. 



(62) 



The nonlinear system in Eq. (62) is solved via Newton iterations. This leads to a linear system 
for each Newton correction, of the form 

J[£n,m+1 — £n,m] = —G{£,n,rn), (63) 

where £„ )m is the mth approximation to £ m . Here J is some approximation to the system's 
Jacobian 

dG dF dF . A ^ 

J = irr^ +a w (64) 

where a = 1/h. In our finite element implementation, we assemble the residual G(£ n ,m) at 
each Newton correction. The implementation of the residual vector is based on the formulation 
presented in Problem 3. However, this formulation makes the determination of the corresponding 
Jacobian rather involved due to the structure of the operator S a ^(w) (see Eq. (42)). Hence, 
we have implemented a Newton-Raphson iteration based on an approximate Jacobian. With 
reference to Theorem 1 and Eq. (43), the Jacobian we assemble is the exact Jacobian of a 
formulation in which the operator product S a ^(w)A 1 (w) is replaced by the operator A a (w, h) 
defined in Eq. (35). In the code accompanying this paper, the final system is solved using the 
direct solver provided by the UMFPACK package (see Davis, 2004). 
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Table 1: Content of the provided zip archive. 



INSTALL . txt: 
Makefile: 
step-f eibm. cc: 
immer sed_f em . prm 
meshes/ : 
prms/ : 
out/ : 
doc/: 



installation instructions; 
standard deal. II makefile; 
main program file; 
default parameter file; 

directory containing a collection of input mesh files, in UCD format; 
examples parameter files; 

empty directory, in which output files will be written; 
Doxygen documentation directory. 



4 Implementation 

4.1 Source files and library requirements 

The included source code is based on the deal .II 7.1.0 library (see Bangerth et al., 2006). In 
what follows, we assume that the user has installed the deal. II library in some directory (we 
have tested our code with deal . II version 7 . 1 and with the latest svn version 7 . 2 pre), and all 
paths will be relative to the install directory path, which we call deal. II. For the program to 
work properly, deal. II should be configured at least with the flag — with-umfpack, to enable 
UMFPACK solver developed in Davis (2004) and used extensively inside the program. 

The provided zip archive should be unzipped under 
deal . II/examples/step-f eibm 

Table 1 provides a summary of the distributed files and directories. If the code is unzipped 
in the above location, it can be compiled by simply typing make at the command line prompt, 
and run with 

./step-f eibm parameter_f ile .prm 

If the file parameter_f ile. prm does not exist, the program creates one with default values, 
which can then be suitably modified by the user should edit for her needs. We distribute all 
the parameter files that were used to produce the results in Section 5 along with the needed 
mesh files. These can be found in the directories prms and meshes, respectively. If the program 
is run without arguments, it is assumed that the problem parameters are those in the file 
parameter_f ile. prm. As mentioned earlier, if the file in question does not already exist, a 
default copy will be created. 

If the user has Doxygen, a complete and browsable documentation of the source code itself 
can be generated in two ways: by typing make in the subdirectory 
deal . II/examples/step-f eibm/doc, 

or by typing make online-doc in the directory deal. II. In the first case, the documentation 
will be accessible at the address 

deal . II/examples/step-f eibm/doc/html/ index .html, 

while in the second case, the command will generate the full deal. II documentation, together 
with the documentation of the step-f eibm source code, which will be accessible at the address 
deal . II/doc/doxygen/deal . II/step_f eibm.html. 

In the second case, the process is a lot longer but the documentation generated thus has the 
added benefits of being fully integrated with the deal. II documentation as well as having 
hyperlinks to all the deal. II classes that have been used in our program. 
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4.2 Parameter and input files 



The behavior of the program is controlled by the ProblemParameters<dim> class, which is 
derived from the deal . II class ParameterHandler and is used to define and to read from a file 
all the problem parameters that the user can set. 

The following is a sample parameter file that can be used with our code. 

# Listing of Parameters 
# 



# Time Stepping 

set Final t =1 

set Delta t = . 1 

set Interval (of time-steps) between output = 1 



# Non linear solver 

set Force J update at step beginning = false 

set Update J cont = false 

set Semi-implicit scheme = true 

set Use spread operator = true 



# Constitutive models available are: INH_0: incompressible Neo-Hookean with 

# P~{e} = mu (F - F~{-T}) ; INH_1: incompressible neo-Hookean with P~{e} = mu F; 

# Circumf erentialFiberModel : incompressible with P~{e} = mu F 

# (e_{\theta} \otimes e_{\theta}) F~{-T>; this is suitable for annular solid 

# comprising inextensible circumferential fibers 



set Solid constitutive model = INH_0 

set Density = 1 

set Viscosity = 1 

set Elastic modulus = 1 

# Dimensional constant for the velocity equation 
set Phi_B = 1 



# Solid mesh information 

set Solid mesh = meshes/solid_square . inp 

set Solid refinement = 1 



# Fluid mesh information 

set Fluid mesh = meshes/f luid_square . inp 

set Fluid refinement = 4 

set All Dirichlet BC = true 

set Dirichlet BC indicator = 1 

set Velocity finite element degree = 2 

# Select between FE_Q (Lagrange finite element space of continuous, piecewise 

# polynomials) or FE_DGP (Discontinuous finite elements based on Legendre 

# polynomials) to approximate the pressure field 
set Finite element for pressure = FE_DGP 
set Fix one dof of p = false 

# Base name used for the output files 

set Output base name = out/square 

# This section is used only when the constitutive model is set to 

# Circumf erentialFiberModel 

subsection Equilibrium Solution of Ring with Circumferential Fibers 
set Any edge length of the (square) control volume = 1. 
set Inner radius of the ring = 0.25 

set Width of the ring = 0.0625 

set x-coordinate of the center of the ring = 0.5 

set y-coordinate of the center of the ring =0.5 

end 
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subsection WO 

# Sometimes it is convenient to use symbolic constants in the expression 

# that describes the function, rather than having to use its numeric value 

# everywhere the constant appears. These values can be defined using this 

# parameter, in the form ' varl=valuel , var2=value2, 
# 

# A typical example would be to set this runtime parameter to 

# 'pi=3. 1415926536' and then use 'pi' in the expression of the actual 

# formula. (That said, for convenience this class actually defines both 'pi' 

# and 'Pi' by default, but you get the idea.) 
set Function constants = 

# The formula that denotes the function you want to evaluate for particular 

# values of the independent variables. This expression may contain any of 

# the usual operations such as addition or multiplication, as well as all of 

# the common functions such as 'sin' or 'cos'. In addition, it may contain 

# expressions like 'if(x>0, 1, -1)' where the expression evaluates to the 

# second argument if the first argument is true, and to the third argument 

# otherwise. For a full overview of possible expressions accepted see the 

# documentation of the f parser library. 
# 

# If the function you are describing represents a vector-valued function 

# with multiple components, then separate the expressions for individual 

# components by a semicolon, 
set Function expression = 0; 

# The name of the variables as they will be used in the function, separated 

# by commas. By default, the names of variables at which the function will 

# be evaluated is 'x' (in Id), 'x,y' (in 2d) or 'x,y,z' (in 3d) for spatial 

# coordinates and 't' for time. You can then use these variable names in 

# your function expression and they will be replaced by the values of these 

# variables at which the function is currently evaluated. However, you can 

# also choose a different set of names for the independent variables at 

# which to evaluate your function expression. For example, if you work in 

# spherical coordinates, you may wish to set this input parameter to 

# 'r ,phi , theta, t ' and then use these variable names in your function 

# expression. 

set Variable names = x,y,t 

end 



subsection force 

set Function constants = 

set Function expression =0; 0; 

set Variable names = x,y,t 

end 



subsection uO 

set Function constants = 

set Function expression =0; 0; 

set Variable names = x,y,t 

end 



subsection ug 

set Function constants = 

set Function expression = if(y>.99, 1, 0); 0; 
set Variable names = x,y,t 

end 
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At the beginning of the parameter we find specifications for the time stepper and for the 
nonlinear solver. In addition, we find information on the constitutive behavior of both the fluid 
and the immersed solid. 

The user can specify the names of the files containing the meshes for the control volume and 
the immersed solid, along with the initial global refinement level for each mesh, in the parameter 
file. In the section pertaining to the control volume, the user can also set the degree of the finite 
element spaces for the fluid velocity as well as the type of the finite element space for the fluid 
pressure. The type and degree of the finite element space for the displacement of the immersed 
domain are automatically set to be the same as those for the velocity of the fluid. A degree 
greater than or equal to two should be selected for the finite element space of the velocity so as 
to ensure proper inf-sup stability. The degree of the pressure space is then automatically set to 
be one less than that for the velocity. 

In the latter part of the parameter file, the user can specify the initial and boundary values 
of the solution as well as the external body forces. Here WO denotes the initial value of the 
displacement of the immersed domain, force denotes the external body force field, uO is the 
initial condition for the velocity and the pressure fields and ug is the Dirichlet boundary condition 
(here configured for a lid-cavity problem). 

The above file, for example, generates the parameters for a lid-cavity problem inside a square 
control volume (read from meshes/f luid_square . inp), with an immersed solid whose mesh is 
given in 

meshes/solid_square . inp. 
4.3 Code structure 

The structure of our program follows closely the structure of most tutorial programs in the 
deal . II library, to which we refer for further explanations and examples. The main class of the 
program is the class ImmersedFEM<dim>, in which all objects and methods to solve the problem 
at hand are defined (including an object of type ProblemParameters<dim>). 

Execution of the solution is triggered in the method run ( ) , which starts the time stepping 
scheme of the DAE system described in Section 3.3, and controls the convergence of the Newton 
iteration scheme for the solution of system Eq. (63). 

Detailed documentation of the code has been embodied in the code itself, and can be auto- 
matically generated with Doxygen. Here we only briefly overview the main ideas behind the use 
of deal. II for immersed methods. 

Due to the nature of the method, two different sets of objects are needed to describe the 
triangulation, the degrees of freedom, etc., of both the fluid and the immersed domains. In 
the code, objects pertaining to the fluid have been denoted with the suffix _f , whereas objects 
pertaining to the immersed solid have been denoted with the suffix _s. For example, tria_s and 
tria_f are the two Triangulation<dim> objects of the solid and fluid domains, respectively. 

In the code, solution vectors and residuals are constructed as 
BlockVector<double> objects and the Jacobian matrix is constructed as a 
BlockSparseMatrix<double> object. This has been done to reflect the logical splitting of these 
entities between the fluid and the solid, and to allow access to the individual blocks at the same 
time. We split the vectors and matrices into two and four parts, respectively. The block vectors 
storing the overall solutions at the current time step and at the previous time step are called xi 
and previous_xi, respectively. The first block of these block vectors pertains to the fluid and 
it is of size n_dofs_up, which is also equal to dh_f .n_dofs(). The second block pertains to the 
solid and has a size of n_dofs_W, which is also equal to dh_s.n_dof s(). 

The various tutorial examples of the deal . II library describe in an exhaustive manner how to 
treat a single triangulation and a single degrees-of-freedom handler for both fluid-only problems 
(e.g., the example program step-35) and elasticity-only problems (e.g., the example program 
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step-44). The most delicate part of immersed methods, however, requires the coupling between 
a fixed background mesh (the fluid), and a moving and deforming foreground mesh (the elastic 
solid) . The deformation of the foreground mesh is achieved very effectively through the 
MappingQEulerian<dim, spacedim> class, which uses the information stored in the displacement 
vector to automatically compute the deformed positions of the mesh and of the quadrature points 
in a Lagrangian way. Notice that while the name suggests an Eulerian description, this object 
in reality performs a Lagrangian iso-parametric transformation from the reference grid, stored 
in tria_s, to the current configuration of the solid via the deformation vector w. Details on 
construction and use of this class are given in Section 4.3.1. 

Evaluation of the quadrature points of the solid on the background fluid mesh is achieved 
through the class FEFieldFunction<dim>, which allows one to evaluate the values of finite 
element fields at arbitrary points. In particular, its method 

FEFieldFunction<dim> : : compute_point_locations is the one that returns the lists required 
to compute the coupling integrals (see Section 3.2) and is used both in the creation of the 
sparsity pattern that features the coupling between the degrees of freedom of the fluid and the 
immersed solid (see Section 4.3.2), as well as in the assembling of the residual vector and the 
Jacobian matrix of the DAE system (see Section 4.3.3). 

4.3.1 Immersed map 

Whenever it is necessary to compute the deformed configuration of the solid, an iso-parametric 
displacement is superimposed on each node of the triangulation of the solid. This process is 
transparent to the user and is performed by the class MappingQEulerian<dim>. In our code, 
we pass an object of this class as an argument to all the standard deal. II classes which are 
involved in computing the finite element values and their gradients on the deformed cells of the 
triangulation of the solid. In the following code snippet we illustrate this process that takes 
place at the beginning of the computation of the residual and of the Jacobian: 

MappingQEulerian<dim> * mapping; 

template <int dim> 
void 

ImmersedFEM<dim> : : residual_and_or_ Jacobian ( . . .) 
{ 

if (mapping != NULL) delete mapping; 

if (par . semi_implicit == true) 

mapping = new MappingQEulerian<dim , Vector<double> , dim> 
(par. degree, previous_xi .block(l) , dh_s) ; 

else 

mapping = new MappingQEulerian<dim, Vector<double> , dim> 
(par. degree, xi.block(l), dh_s) ; 



FEValues<dim,dim> f e_v_s_mapped 



(♦mapping, 
fe_s, 
quad_s, 

update_quadrature_points) ; 



This code snippet illustrates how to instantiate an iso-parametric mapping based on the 
current displacement solution, given by xi.block(l) or on the previous displacement solution 
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previous_xi .block(l) . We refer to the deal. II documentation of the class 
MappingQEulerian for further details on the meaning of each of the arguments passed to the 
constructor of the class. Here it is important to notice that, once a mapping from the reference 
configuration to the deformed configuration is available, it is used in all instantiations of those 
classes which compute the values and the gradients of the basis functions on the deformed 
configuration (i.e., FEValues<dim,dim>). 

Setting the parameter "Semi-implicit scheme" to true in the parameter file (see Sec- 
tion 4.2) will set the variable par . semi _ implicit to true in the above snippet of code. The 
consequence of this choice is that, while the elastic response of the solid is computed at its 
current configuration, i.e., the Piola-Kirchhoff stress is still computed using xi.block(l), the 
body force corresponding to this stress is applied to the fluid surrounding the body at the lo- 
cation previous_xi .block (1) , instead of xi. block (1). In other words, the operator defined 
in Eq. (35), and later split in the change of variable operator and in the Lagrangian elastic 
operator in Theorem 1 (see Eq. (42)), will use xi.block(l) in place of the variable w and 
previous_xi .block(l) in place of the variable h. 

This splitting preserves the consistency of the method, and removes the nonlinearity due to 
the change of variable from the system at the cost of introducing a CFL condition on the time 
stepping scheme (for a more detailed discussion on this topic see Boffi et al., 2007; Heltai, 2008), 
which ceases to be asymptotically stable. 

4.3.2 Sparsity pattern 

A SparsityPattern is a deal. II object which stores the nonzero entries of a sparse matrix. 
Since we are using a BlockSparseMatrix<double> class to store the Jacobian of the DAE 
system, we need a SparsityPattern for each of the sub-blocks of this block. The snippet of 
code that generates the coupling sparsity pattern is given by 

FEFieldFunction<dim, DoFHandler<dim> , Vector<double> > 
up_field (dh_f , tmp_vec_n_dof s_up) ; 

vector< typename DoFHandler<dim> : : active_cell_iterator > cells_f; 

vector< vector< Point< dim > > > qpoints_f; 

vector< vector< unsigned int> > maps; 

vector< unsigned int > dof s_f (f e_f . dof s_per_cell) ; 

vector< unsigned int > dof s_s (f e_s . dof s_per_cell) ; 

typename DoFHandler<dim,dim>: : active_cell_iterator 
cell_s = dh_s .begin_active() , 
endc_s = dh_s.end(); 

FEValues<dim,dim> f e_v_s(immersed_mapping, fe_s, quad_s, 

update_quadrature_points) ; 

CompressedSimpleSparsityPattern spl(n_dof s_up, n_dof s_W) ; 
CompressedSimpleSparsityPattern sp2(n_dofs_W , n_dof s_up) ; 

for(; cell_s != endc_s; ++cell_s) 
{ 

f e_v_s .reinit (cell_s) ; 
cell_s->get_dof .indices (dof s_s) ; 
vector< Point< dim > > &qpoints_s 

= fe_v.get_quadrature_points() ; 

up_f ield. compute_point_locations (qpoints_s , 

cells_f, qpoints_f, maps); 
for(unsigned int c=0; c<cells_f . size() ; ++c) 
{ 
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cells_f [c] ->get_dof _indices (dof s_f ) ; 
for(unsigned int i=0; i<dof s_f . size () ; ++i) 
for (unsigned int j=0; j<dof s_s . size() ; ++ j ) 
{ 

spl .add(dof s_f [i] ,dof s_s [j] ) ; 
sp2.add(dof s_s [j] ,dof s_f [i] ) ; 

} 

} 

} 

Here an FEFieldFunction<dim> object is constructed with a dummy finite element vector 
field (tmp_vec_n_dof s_up) to have access to its member function 

FEFieldFunction<dim> : : compute_point_locations. This member function takes as input 
the location of the quadrature points in each solid cell qpoints_s (computed with the 
FEValues<dim> object f e_v_s, initialized with the mapping described in Section 4.3.1) and fills 
up a series of vectors, which allow the computation of the integrals as explained in Section 3.2. 
These vectors are respectively: 

• cells_f : the vector of all fluid cells containing at least one of the quadrature points of 
the immersed domain; 

• qpoints_f: a vector of the same length as cells_f , containing the custom vector of 
quadrature points in the fluid reference (unit) cell, which gets transformed via the fluid 
mapping to the subset of solid quadrature points qpoints_s (that happen to be in the cell 
in question); 

• maps: a vector of the same length as cells and qpoints_f , which contains vectors of 
indices of the solid quadrature points to which the fluid quadrature points refer to, i.e., 
qpoints_f [i] [j] is mapped by the fluid mapping to the same physical location to which 
the point qpoints_s [maps [i] [j]] is mapped by the solid mapping. 

In the construction of the sparsity patterns, only the first vector, cells_f , is used since we 
only need to know which degrees of freedom are coupled. In particular, all degrees of freedom in 
the fluid cells contained in cells_f will couple with the solid cell identified with the cell iterator 
cell_s. These couplings are computed in the innermost for-loop. 

4.3.3 Residual and Jacobian 

Similarly to what happens for the computation of the sparsity pattern, we use an object of type 
FEFieldFunction<dim> to compute the location of the quadrature points of the immersed solid 
within the fluid cells. Assembly of the coupling matrices is then possible by looping over all 
solid cells, and constructing custom quadrature formulas to use with the fluid cells in order to 
compute the integrals explained in Section 3.2. The following snippet of code explains the most 
relevant points: 

// Loop over solid cells 

for(cell_s = dh_s .begin_active() ; cell_s != endc_s; ++cell_s) 
{ 

fe_v_s_mapped.reinit(cell_s) ; 

up_f ield. compute_point_locations (f e_v_s_mapped.get_quadrature_points() , 

fluid_cells, 
fluid_qpoints, 
f luid_maps) ; 
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// Cycle over all of the fluid cells that happen to contain some of 
// the the quadrature points of the current solid cell, 
for (unsigned int c=0; c<f luid_cells . size () ; ++c) 
{ 

f luid_cells [c] ->get_dof .indices (dof s_f ) ; 

// Local FEValues of the fluid 

Quadrature<dim> local_quad (f luid_qpoints [c] ) ; 

FEValues<dim> local_fe_f_v (fe_f, 

local_quad, 
update_values I 
update_gradients I 
update_hessians) ; 

local_fe_f_v.reinit(fluid_cells[c] ) ; 

// Use the local_fe_f_v as you would normally do: 

for(unsigned int i=0; i<f e_s .dof s_per_cell ; ++i) 
{ 

unsigned int wi = i + f e_f . dof s_per_cell ; 
comp_i = fe_s . system_to_component_index(i) . first ; 
for(unsigned int q=0; q<local_quad. size () ; ++q) 
{ 

unsigned int &qs = f luid_maps [c] [q] ; 



local_res [wi] -= par.Phi_B 

* local_up [q] (comp_i) 

* f e_v_s . shape_value(i ,qs) 

* f e_v_s . JxW(qs) ; 

In the snippet above, we show how the term — f K u(s + w(s,t),t) ■ y(s)ds is assembled in 
practice. The point locations are computed by up_f ield. compute_point_locations. We loop 
over the filled vectors to compute the coupling between each of the fluid cells, f luid_cells [c] , 
and the solid cell cell_s. Since the computed quadrature points in the fluid reference cells 
are not standard (i.e., they are not located at Gauss quadrature points), we need to create a 
custom quadrature formula containing the points of interests (the object local_quad, initialized 
with f luid_qpoints [c] ) as well as an FEValues object, local_f e_f _v, to calculate values and 
gradients of the fluid shape functions at the solid quadrature points. 

These custom FEValues are then initialized with the fluid cell f luid_cells [c] . Notice that 
the correspondence between the indexing in the solid quadrature points and in the fluid custom 
quadrature is given by f luid_maps [c] [q] . The rest follows the standard usage of the deal. II 
library, as can be found in any of the deal. II example programs. 

5 Numerics 

We present in this section two numerical experiments that highlight the aspects of the accuracy 
and error convergence properties as well as the volume conservation feature of our numerical 
method. 

5.1 Static equilibrium of an annular solid comprising circumferential fibers 
and immersed in a stationary fluid 

This numerical test is motivated by the ones presented in Both et al. (2008); Griffith and Luo 
(2012). The objective of this test is to compute the equilibrium state of an initially undeformed 
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thick annular cylinder submerged in a stationary incompressible fluid that is contained in a rigid 
prismatic box having a square cross-section. Our simulation is two-dimensional and comprises 
an annular solid with inner radius R and thickness w, and filled with a stationary fluid that 
is contained in a square box of edge length I (see Fig. 3). In this setting, the reference and 
the deformed configurations of the annular solid can be conveniently described using the polar 
coordinate systems, whose origins coincide with the center of the annulus and whose unit vectors 
are given by (ur,uq) and (u r ,ug), respectively. This ring is located coaxially with respect to 



we, u g 




O 



Figure 3: The reference and deformed configurations of a ring immersed in a square box filled 
with stationary fluid 



that of the box and it is subjected to the hydrostatic pressure of the fluid pi and p a at its inner 
and outer walls, respectively. Negligible body forces act on the system and there is no inflow or 
outflow of fluid across the walls of the box. Since both the solid and the fluid are incompressible, 
it is expected that neither the annulus nor the fluid will move at all. Therefore, the problem 
reduces to determining the equilibrium solution for the Lagrange multiplier field p. The elastic 
behavior of the ring is governed by a continuous distribution of concentric fibers lying in the 
circumferential direction. The first Piola-Kirchhoff stress for the ring is then given by 

P = -p s F~ T + fi e fiie ®ue, (65) 

where /i e is a constant and p s is the Lagrange multiplier that enforces incompressibility of the 
ring. As alluded to earlier, the reference configuration and the deformed configuration of the 
ring must coincide because of incompressibility, and the fact that the deformation of the ring 
must be axisymmetric in nature. For F = I the constitutive response for the Cauchy stress can 
then be written as 

t s = -p s \ + fx e u e ®ue, (66) 

where, for the deformation at hand, ibQ = uq. The balance of linear momentum for the ring can 
be obtained from Eq. (2) as 

-grad(p s ) +n e div(u e ®ug) = 0. (67) 

Noting that grad (p s ) = (dp s /dr)u r + (l/r)(dp s / 'd9)ug, and that div (ug <g> ug) = —(l/r)u r , 
Eq. (67) can be rewritten as 

-7T-^ = ° ^ = 0. (68) 

or r 69 



19 



From Eq. (68), it can be concluded that the Lagrange multiplier enforcing incompressibility p s 
is an axisymmetric function of the form 

p s = c- / u e ln(^), (69) 

where c is a constant. The satisfaction of the traction boundary conditions at the inner and 
outer walls of the ring demand that p s \ r= R = p% and p s \ r =R+w = Po and hence we can obtain 
that 

p s = Pi -fi e ln(^^j , p = Pi - //In (l + *0 (70) 

Note that Lagrange multiplier p defined over the control volume corresponds to p s in the region 
occupied by the solid. By constraining the average value of p over the entire control volume to 
be zero we arrive at the following solution for the equilibrium problem: 



p 



Po = -%({R + wf -R 2 ) for R + w<r, 



p s = ^ e ln(^) - % [(R + w) 2 - R 2 ) for R<r<R + w, (71) 
Pi = n e ln(l + l)-^£ ((R + w) 2 -i? 2 ) for r<R, 



with velocity of fluid u = and the displacement of the solid w = 0. Note that Eq. (71) is 
different from Eq. (69) of Bom et al. (2008), where p varies linearly with r (we believe this to 
be in error). 

For all our numerical simulations we have used R = 0.25 m, w = 0.06250 m, I = 1.0 m and 
\x e = lPa and for these values we obtain pi = 0.16792 Pa and p Q = —0.05522 Pa using Eq. (71). 
We have used p = 1.0kg/m 3 , dynamic viscosity /i = 1.0Pa-s, and time step size h = lxl0~ 3 s 
in our tests. For all our numerical tests we have used Q2 elements to represent w of the solid, 
whereas we have used (i) Q2/P1 elements, and (ii) Q2/Q1 elements to represent v and p over 
the control volume. We present a sample profile of p over the entire control volume and its 
variation along different values of y, after one time step, in Fig. 4 and Fig. 5 for Q2/P1 and 
Q2/Q1 elements, respectively. 

We assess the convergence property of our numerical scheme by obtaining the convergence 
rate of the error between the exact and the numerical solutions of this equilibrium problem. 
The order of the rate of convergence (see, Tables 2 and 3 for Q2/P1 and Q2/Q1 elements, 
respectively) is 2.5 for the L 2 norm of the velocity, 1.5 for the H 1 norm of the velocity and 1.5 
for the L 2 norm of the pressure which matches the rates presented in Boffi et al. (2008). In all 
these numerical tests we have used 1856 cells with 15776 DoFs for the solid. 

The parameter files used for these tests can be found under the directory 
prms/RingEqm_XXX_f ref _Y_param.prm, where XXX is either dgp or feq and Y is 4, 5, 6 or 7, 
according to the type of pressure finite element and to the fluid refinement level. The tests can 
be run under the step-f eibm directory, by typing 
. / step-f eibm prms/RingEqm_XXX_f ref _Y_param . prm 



5.2 Disk entrained in a lid-driven cavity flow 

We test the volume conservation of our numerical method by measuring the change in the area 
of a disk that is entrained in a lid-driven cavity flow of an incompressible, linearly viscous fluid. 
This test is motivated by similar ones presented in Griffith and Luo (2012); Wang and Zhang 
(2010). Referring to Fig. 6, the disk has a radius R = 0.2 m and its center C is initially positioned 
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Figure 4: The values of p after one time step when using PI elements for p 



Table 2: Error convergence rate obtained when using PI element for p after one time step 



No. of cells 


No. of DoFs 


IK - u llo 




IK-u||i 




\\Ph ~p\\o 




256 


2946 


2.00605c-05 




1.95854e-03 




6.71603e-03 




1024 


11522 


3.69389e-06 


2.44 


7.44696e-04 


1.40 


2.47476e-03 


1.44 


4096 


45570 


5.76710e-07 


2.68 


2.25134e-04 


1.73 


8.74728c-04 


1.50 
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Figure 5: The values of p after one time step when using Ql elements for p 
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Table 3: Error convergence rate obtained when using Ql element for p after one time step 



No. of cells No. of DoFs \\uh — u||o || u /i — u l|i \\Ph — p\\o 

256 2467 4.36912e-05 - 2.79237e-03 - 7.393f0e-03 

f024 9539 6.f4959e-06 2.83 9.02397e-04 1.63 2.42394e-03 1.61 

4096 37507 1.28224e-06 2.26 3.49329e-04 1.37 9.10608e-04 1.41 

16384 148739 2.33819e-07 2.46 1.25626e-04 1.48 3.27256c-04 1.48 



at x = 0.6 m and y = 0.5 m in the square cavity whose each edge has the length I = 1.0 m. Body 
forces on the system are negligible. The two different constitutive models for the elastic response 
of the disk which we have used for our simulations are as follows: 

case 1: P = -p s \ + // (F - F~ T ) , (72) 
case 2: P = -p s \ + //F. (73) 

We have used the following parameters: p = 1.0kg/m 3 , dynamic viscosity /i = 0.01 Pa-s, shear 
modulus p e = 0.1 Pa and U = l.Om/s. For our numerical simulations we have used Q2 elements 
to represent w of the disk whereas we have used Q2/P1 element for the fluid. The disk is 
represented using 320 cells with 2626 DoFs and the control volume has 4096 cells and 45570 
DoFs. The time step size h = lxlO -2 s. We consider the time interval < t < 8 s during which 
the disk is lifted from its initial position along the left vertical wall, drawn along underneath the 
lid and finally dragged downwards along the right vertical wall of the cavity (see, Figs. 7 and 
10). As the disk trails beneath the lid, it experiences large shearing deformations (see, Figs. 8 
and 11). Ideally the disk should have retained its original area over the course of time because 
the incompressibility of the media and the nature of the motion require that the disk change its 
shape only and not its volume. However, from our numerical scheme we obtain an area change 
of the disk of about 6% for case 1 (see Fig. 9) and about 4% for case 2 (see Fig. 12). 

The parameter files used for these two tests can be found under the directory 
deal . II/step-f eibm/prms, and are named LDCFlow_Ball_DGP_INHO_param.prm and 
LDCFlow_Ball_DGP_INHl_param.prm respectively. The tests can be run under the 
deal . II/step-f eibm directory, by typing 
./step-f eibm prms/LDCFlow_Ball_DGP_INHO_param.prm 
and 

. /step-f eibm prms/LDCFlow_Ball_DGP_INHl_param.prm 
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Figure 6: The initial configuration of an immersed disk entrained in a flow in a square cavity 
whose lid is driven with a velocity U towards the right 
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Figure 7: The motion of the disk for case 1 at different instants of time 
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Figure 8: Enlarged view of the disk for case 1 depicting its shape and location at various instants 
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Figure 9: The percentage change in the area of the disk for case 1 over time 
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Figure 10: The motion of a disk for case 2 at different instants of time 
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Figure 11: Enlarged view of the disk for case 2 depicting its shape and location at various 
instants of time 
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Figure 12: The percentage change in the area of the disk for case 2 over time 
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